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Abstra ct 

Transonic aeroelastlc solutions based upon 
the transonic small perturbation potential equa- 
tion are studied. Time-fltarching transient solu- 
tions of plunging and pitching airfoils are 
analyzed using a complex exponential modal 
identification technique, and seven alternative 
Integration techniques for the structural equa- 
tions are evaluated. The HYTRAN2 code is used 
to determine transonic flutter boundaries versus 
Mach number and angle-of-attack for NACA 64A0iO 
and MOB A-3 airfoils. In that code, a monotone 
differencing method, which eliminates leading 
edge expansion shocks. Is used to solve the 
potential equation. When the effect of static 
pitching moment upon the angle-of-attack Is 
Included, the MBB A-3 airfoil can have multiple 
flutter speeds at a given Mach number. 
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Introduction 


The calculation of aeroelastlc response 
characteristics In the transonic speed range Is 
of much current interest since the avoidance of 
transonic flutter Is a key design problem. 

Prior to the development of computer programs 
capable of solving transonic aerodynamic pro- 
blems, linear subsonic and supersonic solutions 
were extended Into this regime even though the 
assumptions of the underlying theory were vio- 
lated. The LTRAN21 computer program solves 
the two-dimensional, low frequency, transonic, 
small perturbation potential equation. Its 
application to a simple aeroelastlc stability 
problem was demonstrated by Ballhaus and 
GoorjIan.Z Ref. 2 Illustrates the two ap- 
proaches which have been followed In the appli- 
cation of transonic aerodynamic calculations to 
aeroelasticity ; namely, harmonic analysis and 
time-marching analysis. The former assumes that 
the unsteady aerodynamic forces are locally 
linear and utilizes traditional modal super- 
position of harmonic loads while the latter 
delays the question of linearity by calculating 
the transient response of the coupled aero- 
dynamic-structural system. If the assumption of 
local linearity Is warranted, the harmonic 
analysis approach offers a significant computa- 
tional savings. Rlzzetta^ examined the time- 
marching technique, using a four-point Adams- 
Moulton Integrator for the structural equations 
of motion and the LTRAN2 code to calculate the 
unsteady airloads. The Initial conditions 
chosen were large enough that significant 
nonlinear effscts occurred In the calculated 
unsteady airloads, Indicating large shock 
motions. Yang, Guruswamy, and Strlz^ studied 
a tlme^narchlng scheme which coupled transient 
LTRAN2 solutions and a structural Integrator 
which assumed a linear variation of velocities. 
They also compared harmonic flutter analyses^ 
using oscillatory airloads derived from LTRAN2 
and UTRANS26. 
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■itn^v' tHo ui i’jinal I TRAN,' i WiJi' aiiuf’ati* only 
at low frotiuonolos (K<,0?f)), several tntprove- 
ments were Mile to inoreOSC its ranye of appll- 
tabllUy, ilouwlnk and van der Vooren' studied 
the efteot of retainintj the time derivative 
terms in the boundary and auxiliary conditions 
in their iTRANi’.NlR rode and claimed accuracy to 
k . n.a. Rfi?ettd and ChinH retained, in 
addition, the h1>jh frequency (ttt) 1’’ 
yovernlni) equation. Isoyal^ has included both 
of these effects In a computer code which has 
been used to generate oscillatory airloads for a 
harmonic flutter analysis of a two-dimensional 
airfoil section. The analysis demonstrated a 
significant transonic dip In the section's 
flutter boundary. Borland and RlzzuttalO have 
developed a three-dimensional unsteady aero- 
dynamic code, XTRAN3S, and have utilized a cen- 
tered difference structural Integration tB-onl- 
que to obtain transient time -marching flutter 
solutions. 

Application of these time dependent 
small -perturbation codes has been hampered by 
numerical stability problems which occur for 
large amplitude motion and/or large 
angles-of-attack, The current study modifies 
the LTRAN2-NLR Code with the monotone 
differencing scheme of Engquist and Osherli 
In order to eliminate the source of this 
problem, namely expansion shocks near the 
airfoil leading-edge. The resulting code Is 
termed HTTRAN2. This code has been used to 
study the accuracy and stability o,' the various 
numeHeaT Integration techniques which may be 
used In transient time-marching calculations. 

The accuracy of the finite-difference computer 
program algorithms may be established using 
recent analytic solutions by Bland^^ of the 
linearized tTRANZ, HYTRAN2, and EXTRAN2 
aerodynamic problems. The linearized version 
of the HYTRAN2 code Is then used to Investigate 
the accuracy of various numerical Integration 
techniques for the structural equations of 
motion. A comparison Is made of transonic 
flutter boundaries of the NACA 64A010 airfoil 
calculated by four different nonlinear 
transonic codes. The effect of angle-of-attack 
upon the flutter boundaries of the NACA 64A010 
and the MBB A-3 airfoils Is also danonstrated. 
Finally, the effect upon the flutter boundary of 
Including the aoroelastlc twisting resulting 
from the steady pitching moment Is demonstrated. 


U nsteady Transonic Snail Perturbation Equation 

The calculations described herein were ob- 
tained from a modified version of the LTRAN2-NLR 
COde^ which solves the low frequency transonic 
small perturbation (TSP) potential equation 

[1-M^{y*+1)n 2^J ■fxx ^zz ■ 

subject to the following airfoil boundary and 
wake conditions 

tj, * f^ fj. 1 z « + 0. 0 < X < 1 (2a) 

[♦xi *01^” (2b) 


Aerodynamic loads are given in terms of a 
nondlmenslonal pressure coefficient, r.,, 
which in the small perturbation limit becomes 

Up ^ 4'j.) (2c) 


The airfoil coordinates are given by 
z = f(x,t), and T* • 2-{2 -y)M« 2. The 
original ITRANZ' grid of 99 by 79 points 1n 
the X and z direction, respectively, with 33 
points on the airfoil chord was used. The 
original LTRAN2-NLR and LTRAN2 codes are very 
sensitive to angle-of-attack changes and large 
airfoil pitching or plunging motions, and 
numerical Instabilities generated at the loading 
edgel3 can lead to program failure. The 
monotone differencing method of Engquist and 
Osher'l eliminates the leading edge expansion 
shocks which cause this problem and has been 
Incorporated Into the code used for this study 
which Is tenned IIYTRAN2. Solutions obtained 
with the time derivative terms of (2) deleted 
are referred to as LJRAN2 solutions while 
solutions with a term added to (1) 

are referred to as EXTRAN2 solutions. Refs. 8 
and 9 give results obtained from EXTRAN2 type 
codes. 

The frequency limits of LTRAN2 and HYTRAN2 
codes are generally accepted as k-0,075 and 
k**0.4 respectively. These estimates are based 
upon comparisons of results from linearized 
versions of these finite-difference codes 
(obtained by deleting the ij>x'i'xx term In (1)) 
with classical solutions of the subsonic wave 
equation. This has left open the questions of 
convergence of solutions for a given grid and 
accuracy of the solutions. BlandlZ has 
modified the kernel function of the Posslo 
Integral equation to enable computation of exact 
analytic solutions of the LTRAN2 and HYTRAN2 
linearized problems. Figure 1 gives comparisons 
of his exact c» solutions at M 0.8 for 
the three different kernels with LTRAN2 and 
HYTRAN2 results for reduced frequencies up to 
0,5, The finite-difference results were 
obtained using 3 cycles of time marching 
calculations with 360 steps per cycle. The 
agreement between the finite-difference and 
analytical solutions Indicates that the 
finite-difference code Is quite accurate for the 
grid and step size used. Also, the departure of 
LTRAN2 from the analytic EXTRAN2 results above 
k»0.05 confirms the limitations of LTRAN2 
mentioned above. The HYTRAN2 results are In 
better agreement with EXTRAN2 and show a gradual 
departure from EXTRAN2. The selection of an 
upper frequency limit on HYTRAN2 is somewhat 
arbitrary. 


T1me-March1ng Aeroelastic Solutions 

The classical description of a two- 
dimensional, pitching and plunging airfoil Is 
assumed. The airfoil lies between _+l on the 

x-axIs; plunge, h, and lift coefficient, Cj 
are measured positive downwards at the elastic 
axis, a; and pitch, «, and moment coefficient 
about a, C|n, are positive nose up. The equa- 
tions of motion are written In vector notation 
as 
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The static load coefficients, cj. and 
u,« , are subtracted from the total coaffU 
cients during the tlwe-marchlng calculation. 

Thus h and « represent perturbations about 
assumed static operating conditions, The static 
angle-of*ettack is a separate Input to HYTRAN2 
and contributes to fx In (2a), The speed 
Index V « l)/bw,/u determines the density 
(altitude) assumed for a calculation. As V 
Increases fro<n lero, transient responses are 
Initially damped. Further Increases in V will 
usually lead to neutrally stable oscillations 
characterlailng a flutter boundary, Tlte value of 
V at flutter Is termed the flutter speed Index, 
Vf. A fourth order linear state equation may 
be developed from (3) as 

X ■ M + Hu (4) 


(4) while the second term Is a convolution 
Integral giving the forced response. For use as 
a structural integrator In aeroelastlc time 
marching solutions, (5) Is rewritten to reflect 
the evolution of the structural state from cnee 
step n to tltee step nfl, 

X [(n*-l)T] . a(T)x(nT) 


(n+l)T 

+ / exp[A((n+l)T-t)]Bu(t) dt (6) 
nT 


where T Is the step site, 

The alternating-direction implicit solution 
algorlthn used In HYTRANa requires tnrec pieces 
of data to calculate the flow field potential, 
^n+li «t time step n+li l)+n, the potential 
at time n, g)xni the boundary condition at 
time n, and 3)X|^+i, the boundary condition at 
time n+1. Since the lift and manent are not 
known over the Interval nT<t<(n+l)T, the 
Integral In (6) must be approximated, The simp- 
lest approximation Is to assume that 
u(t) « u(nT) over the Interval. Then llu(nT) may 
be taken out of the Integral ami a slight change 
of notation gives 

4+1 • % »fi''n (7) 

where the Integral of the transltlof) matrix Is 


where 

J . ly' yT] . 

|h a ti 
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Structural Integration Technique 

published studies of tlmo-tnarchlng aero- 
el astlc solutions have coupled the structural 
equations of motion to aerodynamic codes using a 
variety of numerical Integration techniques, 
R1zictta3 used as, Adams-ktoulton predictor- 
corrector scheme, Ref. 4 reports an algorltlint 
based upon an assumed linear variation of accel- 
eration while a contered-difference Integration 
technique was used in ref. 10, None of these 
approaches lakes advantage of the linear 
structure of (4), 

Since (4) Is a finite-dimensional linear 
differential equation, its solutlonl*' Is 

t 

x(t) » 4'(t)x(0) +/ exp[A(t-t) ]Bu(t) dt (5) 
” 0 

The state transition matrix 4>(t) o oxp[At] nwy 
be calculated to any desired accuracy as the sum 
of the first n terms of the series expansion of 
the matrix exponential function. The first term 
In (5) Is the homogeneous response portion of 


0 *< / exp[A(T-T)] dt (8) 

0 

An Improvement upon this approximation may 
bo obtained by considering u to vary linearly 
between U|i and u„+i In (6), Then, for small 
time steps, T, the Integ’al Is nearly equal to 
OB(u.,+i + Un)/2, However, u,,^.! Is not 
available at this step of the algorithm and an 
estimate of Un+i * Un + (un - un.i) will 
be used. Ttie resulting algorithm Is 

4+1 • '*'4 ^ Q‘*(3“n““n-l)^2 (9) 

The Integration iiatrkes 4> and o were calculated 
using the program described In ref. 15. Unless 
othttrwlso noted, the results presented wore ob- 
tained froin transient responses calculated using 
(9) for the structural Integrator. After the 
steady state flow field for the static 
angle-of-attack was Obtained the transient was 
excited by a One percent chord displacement 
initial condition on the plunge coordinate. The 
transient was calculated for H50 time steps 
which yielded 3-6 cycles of oscillation of the 
dominant flutter mode for the examples studied. 
To determine a flutter point, several transients 
were calculated for a range of speed Indices, 

V, Typically, speed Indices were used which 
gave slightly subcrltlcal d.amped response and 
slIgK.ty supercritical diverging response. The 
flutter speed Index was then determined by 
interpolation, Once the general nature of the 
flutter boundary was understood, additional 
flutter points could generally be obtained 1n 
this manner with the calculation of two 
responses per flutter point. 
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HHtiijH 5i tural Paraflieters 

Tub structural paranteters for the two 
t'xaaple cases studied are given in Table 1, 
txample 1 Is the as Case A of isoyais*'^ 
which was salectud to have plunging and pitch- 
ing normal modes similar to those of a stream- 
wise section near the wlngtlp of a sweptback 
wl'.'j. The pivotal point for the lowest fre- 
gueniy plunging mode Is 1.44 chord lengths 
ahead of the leading -edge and for the higher 
frequency pitch mode It Is 0. 068 chordlengths 
ahead of midchord. Isogallo has shown a 
significant transonic dip In the flutter speed 
Index of this section with the plunging mode 
becoming the flutter ntode at frequencies ranging 
from 80-1 BO rad/soG, 

Example 2 was chosen to Investigate the 
situation In which the normal mode frequencies 
were close together and for which the higher 
frequency pitch mode would become the flutter 
mode. The small negative static unbalance 
leads to very little Inertial coupling of the 
modes and also causes the section to be 
nearly statically unstable. 

Table 1. SECTION STRUCTURAL PARAMETERS 



Example 1 

Example 2 

a 

-2 

-0.042 

*0 

1.8 

-0.036 

r„2 

3.48 

1.872 

u 

60 

60 

wn. r/s 

100 

23.5 

r/S 

100 

35 

wind off plunge freq., r/s 

71.34 

23.497 

wind off pitch freq., r/s 

533.0 

35.037 


Modal Identification Technique 

To efficiently use the Information con- 
tained in the transients, the least squares 
curve fitting program of Bennett and 
Oesmaralsl' was used. Both h and a responses 
were fit Independently by the function 


m -(jjt 

“ a^j + I e a^cos(ujjt + \(i^) (10) 

1^1 

For the two degree-of-freedom examples treated 
m was set equal to two. The complex modes thus 
obtained, si » oi + jmj, are estimates of 
the eigenvalues of the aeroelastlc system and 
can be plotted as a function of speed Index, 

V, giving s-plane root 1ocl of the coupled 
plunge and pitch modes. Figure 2 gives a 
typical transient response of the MSB A-3 airfoil 
for Example 1 at M » 0.8 and an arvgle-of-attack 


.« j 0. Shown at the top of figure are the h 
and a responses and the fits of the respons'es. 
At the bottom of the figure the time histories 
of the component modes of the fits are shown, 
The fits are very good and Indicate that non- 
linear effects are very small. 

Table 2 g^ves the estimates of the two 
modes for each of two response time hlstorlts. 
Standard deviations of the estimates are com- 
puted from the residuals and are given In 
brackets beneath the estimates. The estimates 
given by the Independent fits are very consis- 
tent with each other although the standard 
deviations are somewhat optimistic. These 
results Indicate that locally linear techniques 
are applicable for this case Involving plunging 
oscillations of ^0.008 semichords and 
angle-of-attack oscillations of +0.25 deg. 


TABLE 2. MODAL ESTIMATES FOR MSB A-3 RESPONSE 
AT M <• 0.8 AND « « 0 DEG. 



si, rad/scc 

$ 2 , rad/sec 

h-fit 

-.12+186.91 

-17.24+J535.S4 


(.OS+j .05) 

{.31+j .32) 

o-fit 

-.02+J86.74 

-17.94+jS36.8l 


(.07+J .07) 

(,42+J .43) 


Airfoils 

Figure 3 shows the profiles of the two 
airfoil studied. Coordinates for the NACA 
64A010 symmetrical airfoil were taken froai 
Abbott and von Doenhoffl^ and for the MBB A-3 
supercritical airfoil from Bland. 19 
Applications of TSP cod-'s to both airfoils have 
been extensively reported (e.g. refs. 3, 8, 20, 
21 for the NACA 64A010 and refs. 21, 22 for the 
MBB A-3), The MBB A-3 theoretical design 
condition Is M » 0.76S, o « 1.3 deg., and 
cj *« 0.S8 and experimental pressure 
distributions are given by Bucclantlnl et 
al23. The airfoil slopes required by HYTRAN2 
were generated using the geometry processor of 
LTRAN2.1 

The MBB A-3 slopes used In ref, 21 are from 
ref. 24 and were obtained from least-squares 
polynomial curve fits to airfoil ordinates. The 
resulting slopes are Inaccurate, particularly 
for the leading-edge region of the upper surface 
and the trail ing-edge region of the lower sur- 
face where the calculated pressures” are 
significantly different from experiment and de- 
sign. 23 This problem has been corrected in 
ref. 22, 


Linearized Aeroelastlc Model 

Pitching and plunging airload coefficients 
from Bland's exact solutions of the linearized 
LTRAN2, HYTRAN2, and EXTRAN2 problems were used 
to generate matrix Pade approxlmants of the 
loads using the technique described by 
Edwards. 25 The approxlmants were then avail- 
able as linear differential equations which were 


4 



coupled to (4) to allow llneap eigenvalue 
anel;^sU of stability as a function of speed 
Index. The approxlmants give a good 
approximation of the loads along the s>ptane 
Imaginary axis and, by analytic continuation, 
are also valid for complex values of s near the 
Imaginary axis. These models, denoted PUTRAN2, 
PIITTRANZ, and PEXTRAN2, were used to evaluate 
time nvirching solutions from the linearized TSP 
equations for sub.:r 1 t 1 cal , critical and 
supercritical flutter conditions. 


Results and Discussion 

Comparison of Time-Marching and Linearized 
Wutlpt^ 

A eroelasttc root loci # Figure 4 compares 
the pTunge ^d pitch mode root locations 
obtained from curve fits of the plunge response 
with root loci calculated using the linearized 
rade models. The time-marching calculations 
were obtained from the linearized LTRAN2 and 
HTRAN2 codes using a flat plate airfoil. The 
Mach number 1$ 0.8, and the Example 1 structural 
model was used with a time step of 0.001 sec. 
Time-marching results are shotm for V « 0.87, 
1.37, and 1.7S corresponding to flutter points 
of PLTRAN2, PEXTRAN2, and an unstable 
condition. These speed conditions for the Fade 
models are Indicated by tick marks. The true 
flutter point for this case Is given by PEXYRAN2 
as Vf » 1.37, uj! " 135 rad/sec, and kf • 

0.127 which Is closely approximated by both the 
PIIYTRANZ and IIYTRAN2 results. While the flutter 
point predicted for this case by PLTRANZ at Vf 
» 0.87, tiif » 101 rad/sec, and kf « 0.15 Is 

reasonably reproduced by the ITRAN2 result, both 
are over 35< low In Vf and significantly low 
In frequency. There is also good agreement 
between the time-marching results and the pade 
model results for both the stable and unstable 
conditions. The pitch mode results do not show 
as good correspondence. This Is due to: 1)The 
matrix Pade approxlmants arc less accurate at 
the high reduced frequencies and damping ratios 
Involved, and 2)t1me-march1ng Integration errors 
are of concern for the relatively high pitch 
frequency of this example, Still, the HYTRAN2 
results are 1 n reasonable agreement with the 
PEXTRAN2 results. These results give confidence 
in the use of time-qnarching aeroelastic studies 
based upon the HYTRAN2 code. 

Structural Equation Integrator Stu^ - 
Table 3 lists seven Integration teclinT^ues which 
may be used for the structural equations of 
motion. Each Integrator was evaluated for Its 
ability to accurately reproduce the PHYTRAN2 
flutter point of fig. 4 as f function of time 
step size, integrators II and 12 are the state 
transition matrix Integrators of (7) and (9) 
while 13 Is (g) Implemented as a predictor- 
corrector. Integrator 14 is a centered 
difference (In time) algoritlw based upon (4) 
while 15 Is the centered difference algorithm 
used by Borland and Rizzettalb based upon 
(3). Integrator 16 has been used by Yang 
et al 21 » 2 Zwhiie 17 1 $ the same algorithm 
Implemented as predictor-corrector and which was 
used In the XTRAN3S code. 26 


Table 4 gives the flutter (plunge) and 
pitch modes Identified from the plunge transient 
responses of Example 1 at Mach 0,8 for s'lep 
sizes of 0.0001-0,003 seconds and a speed Index 
of V <• 1.5. Table 5 gives the two modes 
Identified from the plunge transient responuns 
of Example 2 at M»ch 0.8 for step sizes of 
0.001*0.03 seconds and a speed Index of V » 
0,747, In this case the flutter root evolves 
from the higher frequency pitch mode. Also 
given are eigenvalues of the respective cases 
from PIIYTRANZ which are the "exact" values for 
the two flutter roots. This was ensured by 
forcing the matrix Pade approxlmants to match 
81and'$12 exact load values at the respective 
flutter frequencies. The flutter reduced 
frequencies for the two examples of Tables 4 and 
5 are kf * 0.126 and 0.146, respectively. 

TABLE 3. STRUCTURAL EQUATION INTEGRATORS 
M - 4X„ + OBu„ 

12 ^+1 " + .50B[3u„ - u„^,] 

13 predictor; - #x^ +.5oB[3u„ - u„.|] 

corrector; x^,^^ - ax^ +.5oB[u„^,+ uj 

H Xfl,., - Xfl.i + 2T[Ax^ + BuJ 

16 y,^, - [21 - T^M'^Kly^ -y„., + T^M'^u^ 

16 - lMfT2K/6j‘’[-Kyj,-TKy^-T2l^'j,/3 + uJ 

+ yj/2 

»n .1 - 'n * 

17 predictor: Integrator 16 

corrector: Integrator 16 with u^ replaced 

%4-l 

All of the Integrators, with the exception 
of 16 and II, were accurate for both examples 
for the smallest time steps. While the modes 
Identified by II and 16 for the smallest time 
steps are perhaps acceptable, both of these 
integrators degrade quickly with large time 
steps. Neither predictor-corrector, 13 and 17 , 
was stable for larger time step sizes than the 
Integrator upon which they were based (12 and 
16) although 17 did Improve the accuracy of the 
flutter root estimate over 16. The flutter nmde 
damping estimates of IS, 16, and 17 show a dis- 
concerting behavior of reversing trends as the 
time step Is Increased, whereas the Integrators 
based upon the transition matrix (U, 12 , and 
13) give uniform damping estimate degradation. 

The heavily damped pitch mode of Example 1 Is In 
close agreement with the f'HYTRANZ value for 
small time steps, whereas the plunge mode damp- 
ing value for Example 2 which Is given by 
PHYTRAN2 is only one-half the value estimated 
from the HYTRAN2 results. As the time step 
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TABU 4, AtROIlASTU; HODlb lOtNOUtD TOK feXAWU I at V • 1 .5 and H • 0,8 


a), flutter (plunge) node 






f.sec, 



integrator 

” PHYTRAN2 

.0001 

,0005 

.001 

.002 

",ooT 


' Voi+Tl 46,27 ‘ 

,60Vjl46,59 

2, 63+3 145.77 

4,76+3145.35 

8.68+3144.18 

i2.lT+3r4?3 

12 

,12+3146.66 

,31+3140,20 

,23+3146.54 

,33+3147.88 

.91+3149.91 

13 


,12+Jl 46.65 

.28+3146.03 

unstable 


14 


.13+3146,66 

unstable 




iu 


.13+3146.65 

.26+3146.1 3 

.00+3146.30 

-.69+3147.16 

unstable 

16 


1.07+3146.52 

4.85+3145.37 

8.89+3144.30 

unstable 


17 


.11+3146.66 

,21+3146.07 

-.22+3146.01 

unstable 





b). pitch inode 


'm— j.-*£ a:::w-aTS' » » 

11 

TT'Tt 3»‘ 

-101.32+3498,02 

-97,77+3509.73 

-108.94+3521 .06 

-1 23.91+3534.60 

-146.32+3570.70 

-154.30+3615.40 

12 


-94.78+3507.09 

-100.60+3506.28 

-113,44+3502.19 

-163.45+3503.05 

-283,03+3578.20 

13 


-94.65+3507.16 

-97.10+3507.86 

unstable 



14 


-94.57+3507.61 

unstable 




15 


-94.50+3507.43 

-98,43+3509.76 

-107.65+3512,72 

-1 56.20+3528,83 

unstable 

16 


-99.20+3512.60 

-119.17+3547.47 

-74.33+3614.32 

unstable 


17 


-95.04+3507.37 

-98.57+3517.91 

-73.18+3549.55 

unstable 



TABLE 5. AEROELASTIC MODES IDENTIFIED FOR EXAMPLE 2 AT V - .747 and M « 0.8. 


a), flutter (pitch) mode 





T, sec. 

Integrator 

PHYTRAN2 

.001 

.005 

.01 

.02 

.03 

11 

.02+329.57 

-.02+329.37 

-.12+329.23* 

-.24+328.99 

-.38+328.37 

-.40+327.78 

12 


-.00+329.43 

.00+329.49 

.06+329.48 

,19+329.22 

.28+328.83 

13 


-.01+329.43 

-.03+329.51 

-.08+329.58 

unstable 


14 


-.00+329.42 

.08+329.89 

unstable 



15 


-.00+329.42 

-.01+329.54 

,03+329.70 

.16+330.08 

unstable 

16 


-.03+329.31 

-.20+329.14 

-.31+329.12 

.89+329.45 

unstable 

17 


-.00+329.42 

-.03+329.73 

.098+330.55 

unstable 



b). plunge mode 




11 

-7.46+325,02 

-15.73+322.31 

-15.94+322.53 

-16.40+322.83 

-20.20+322.40 

-28.67+318.10 

12 


-15.40+322.31 

-14.49+322.34 

-13.97+321.96 

-15.00+320.77 

-17,80+319.06 

13 


-15.40+322.32 

-14.36+322.46 

-13.34+322.45 

unstable 


14 


-15.46+322.35 

-11.10+322.84 

unstable 



15 


-15.47+322.34 

-14,58+322.65 

-14.04+322.68 

-14.38+322,64 

unstable 

16 


-16.05+322.21 

-17.46+322.37 

-19.50+323.58 

- 9.32+336.49 

unstable 

17 

1 , ,, 1 

-15.40+322.24 

-14.1 7+322.36 

-12.59+322.78 

unstable 
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was increased, all of the Integrators eventual ly 
experienced a high frequency Instability which 
could not be attributed to a structural mode 
instability. Integrators 11 and 12 were 
unstable for tiwe steps of .004 sec. for 
Example 1 and .04 sec. for Example 2, At these 
time steps, the structural Integration is being 
performed with only 3-5 steps per cycle of 
oscillation of the higher frequency mode. The 
Integrator giving the best overall performance 
was 12 which gave usable flutter mode damping 
results for time steps 50t larger than the noxt 
best Integrator, 15. Integrator 12 was used for 
the remainder of the study. 

Transo ni c Flutter Anal ysis 

111 this section, the techniques demonstra- 
ted above are applied to the full nonlinear 
equation, (1). Flutter boundaries are presented 
for the structural dynamic nvidel of Example 1 of 
Table I and for NACA 64A0lu ind MBB A-3 air- 
foils. Steady-pressure distributions for the 
ranges of Hach number and angle-of-attack inves- 
tigated are given in figures 5 and 6, In the 
cases where the steady shocks are located near 
the trailing edge, the shock strengths are 
of concern, For these cases In figures 5 and 6, 
the computed Mach numbers ahead of the upper 
surface shocks are only slightly less than 1.3 « 
the Mach number at which entropy rises start to 
become si gnf leant Indicating that these 
calculations are on the edge of and possibly 
outside the range of applicability of TSP 
theory. At a .■ -1,5°, the Mach number on the 
lower surface of the MBB A-3 is such that these 
cases may be beyond the scope of the TSP theory. 


Effect of Integration Step Size - The 
effect of i ntegration step si ze upon the modes 
identified from the angle-of-attack response is 
given in Table 6 for the NACA 64A010 at 
F' ■ 0.8. The flutter speed index is Vf ■ 1.07 
anv kf • 0.139. A comparison with the 
Intv'.rator 12 results of Table 4 indicates 
sligu'.iy smaller changes due to step size in the 
mode e.,,imates for the nonlinear case and gives 
confider s in the application of the 
time-marching technique to nonlinear aeroelastic 
problems. 


Table 6, - MODAL ESTIMATES FOR NACA 64A010 AT 
M » 0.8, n - 0 DEG,, AND Vf ■ 1.07 AS A 
FUNCTION OF T. 



T, sec 

SI. r/s 

sg, r/s 


.0001 

.11+J115.21 

-31.89+3536.26 


.0005 

.08+jl 15.09 

-32.95+3536.17 


.001 

.02+J115.18 

-35.76+3537.04 


.002 

-.08+3115.54 

-42.03+3544.83 


.003 

-.07+3116,08 

-46.07+3551.90 


Effect of Anplitude - The effect of ampli- 
tude on We resyianse H flutter was studied for 
the MBB A-3 alrfoli three conditions! 


l)M "O.TTS and n » 0 deg. where the flow is 
subcritical, 2)M * 0.7937S and a » 0.5 deg. 
where there is a strong shock aft of midchord, 
and 3)M » 0.8 and n » 0.5 deg. where the Shock 
Is at the traiUng-edge. Table 7 gives the 
modal estimates including the amplitude 
estimates from the «»fit for M ■ 0.775 which is 
typical. The flutter speed index is Vf ■ 1,21 
and kf » 0.125. There is a surprisingly small 
effect of amplitude upon the damping of the 
flutter mode. The damping ratio, c, changes 
from +0.011 for ho • 0.01 where the 
oscillation amplitude is +0.28 deg, to -0,004 
for ho • 0.10 where the dscillation amplitude 
is +2.97 deg. Inspection of the pressure 
distribution for the ho « 0.1 case during a 
cycle of oscillation revealed a complex pattern 
of shock motion. During a portion of the cycle, 
a strong shock forms near the 70jt chord 
position, travels forward, weakens, and 
disappears between 40-50iI chord, This is the 
type B shock motion studied by Tidjeman,27 
Also, duping the nose-down portion of the cycle 
a strong shock forms ut 5-1 0> chord on the lower 
surface, It is interesting that the modal 
amplitudes, ai, are nearly linear with respect 
to ho, even for the large amplitudes studied. 

At M « 0.79375 the effect Of amplitude upon the 
flutter mode was smaller than that shown In 
Table 7, the change in the flutter root in going 
from ho » 0.01 to 0.10 being s •> -0.30+J0.69 
rad/sec. At M ■ 0.8 the corresponding change 
was s ■ -1.41 -JO. 19, In all cases the standard 
deviations of the estimates, which did not vary 
appreciably with hn, were small and of the 
order given in Table 1 and the amplitudes varied 
nearly linearly with ho, Calculations of 
amplitude effects were made for the NACA 64A010 
and similar small effects on the modal estimates 
were observed. These results contrast with 
Dowell et al28 who studied forced oscillations 
of the NACA 64A006 airfoil using LTRAN2 and 
postulated that nonlinear aerodynamics would be 
Important at k 0.1 for oscillatory amplitudes 
greater than 0.5 degrees, In summary, it 
appears that, for a given Mach number and steady 
angle-of-attack, the aeroelastic response of 
airfoils with dynamics similar to those studied 
may be treated as locally linear in amplitude 
within the limits of small disturbance theory. 

TABLE 7. MODAL ESTIMATES FOR MBB A-3 AT 
M • 0.775 AND 0 ■ 0 DEG. AS A FUNCTION OF ho 


h 

0 

s . r/s 
1 

s , r/s 
2 

a ,deg 
1 

a ,deg 
2 

.01 

.02 

.04 

.10 

-1.31+3117.17 
-1.26+3117.1 7 
-1. 02+3117.15 
.42+3117.24 

-46.53+3529.55 

-46,20+3528.77 

-46.48+3528.28 

-48.85+3527.14 

.28 

.60 

1.21 

2.97 

.34 

.66 

1.31 

3.24 


Comparison of TSP Flutter Boundaries For 
Example 1 - The transonic fltjtter boundary of 
Example I with an NACA 64A010 airfoil at 
a w 0 deg. has been studied by Isogal9.18 
using an EXTRAN2 code to obtain harmonic 
perturbation airloads and by Ehlers and 
Weather1ll29 using a transonic code, 0PTRAN2, 
to compute the hannonic linear perturbation 
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fluth'p in luiUi 

worn then lalwulatol usiiuj im.nttuniil 
frniiii('my.=act«aln techniques. Ihe tomparlsoii of 
these flutter hound, I'-ies with that obtained from 
t1t»ew»8Pchlng solutions using HYTRAN2 is shown 
ill figure 7. Also shown are two flutter points 
at Math O.ft and 0.82S which used LTRAN? 
aerodynamics. As in figure 4, the ITRAN? 
flutter I’olnt is conservative by up to 60% in 
flutter speed imiex while the IIYTRAN2 results 
are slightly unconservative below M » 0.85 when 
et^mpared to EXTRAN2. The HYTRAN2 results 
compare favorably with Isogai's EKTRAN2 results 
with both giving minimun flutter speed indices 
of Vf ^ 0.5, The HYTRAN2 and EXTRAN2 curves 
are displaced from each other by M » 0,015 which 
may be due to differences In steady pressure 
distribution caused by grid differences. The 
minimum Vf obtained using HYTRAN2 occurs at M 
« 0.85 where the steady shock 1s at x/2b - 
0.75. The flutter speed Index Is slightly 
greater at M » 0*875 where isogal was unable to 
obtain a flutter solution, The 0PTRAN2 results 
(using NACA 64A010A19 airfoil coordinates) are 
In reasonable agreement with HYTRAN2 and EXTRAN2 
for M < 0,82. Multiple flutter points are 
predicted by 0PTRAN2 for Mach nunbers above the 
minimum Mf (0.85 < M < 0.87) and the HVT0AN2 
calculations have confirmed this effect at M • 
0.875. Ref, 29 indicates a complex flutter 
boundary at larger values of V for 
0.88 < M <0.90 which has not been studied with 
HYTRAN2, Thus the flutter boundary between M « 
0.875 and M >• 0,9 1s not shown, For M > 0.9 
the shock has moved off the trail Ing-ecfge and 
the EXTRAN2 and HYTRANg results are In good 
agreement. 

Effect of Anal e-of-Attack « The angle^of- 
attacFTs**Kowh to WTh^Fportant parameter 
affecting transonic flutter. AshleyoO lists 
several Instances of such an effect and Doggett 
and Rlcketts^l have studied the effect of 
angle>of-attack urron an arrow-wing configura- 
tion. Edwards32 gives subcrltlcal damping es- 
timates Indicating changes In damping ratio of 
0,02 for a 0.3 deg Increase In angle-of-attack 
for a supercritical wing. Yates et al33 give 
wind tunnel flutter test results of a similar 
wing at several different angles-Of-attack, 

A novel feature of the flutter boundaries shown 
Is the occurence of multiple flutter speed 
Indices for Mach numbers less than that at the 
bottom of the transonic dip. That is, as the 
flutter speed Index decreased with Increasing 
Mach number, cases were studied In which further 
decreases In speed Index resulted in the flutter 
Mach number decreasing also, This curl -back of 
the flutter boundary occured for angle-of-attack 
changes of approximately two degrees. Houwlnk 
et al34 report a similar occurrence. 


more gradual Steepening than those of fig. 9 as 
Mach number increases. Comparing figs. 6 and 9 
indicates that the abrupt steepening of the MSB 
A-3 boundaries occurs at the Mach number at 
which the upper surface shock forms . The 
minimum Vf of both airfoils is approximately 
0.5 and is not a strong function of a. Also, 
the width of the transonic dip near the minimum 
Vf Is greater for the NACA 64A010 than for the 
M8B A-3, which correlates with the change in 
Mach number required for the shock to travel 
from near midchord to the trailing edge (see 
figs. 5 and 5), For both airfoils, the boundary 
rises sharply when the shock reaches the 
trailing edge, 

Figure 10 gives the bending mode root locus 
for the NACA 64A010 versus V for several angles- 
of-attack and Mach numbers. Fig. 10a presents 
root loci typical of the case In which the shock 
has not yet mov,',d aft on the airfoil (see fig, 
5). Increasing angle-of-attack causes a loss of 
damping and a drop In flutter speed index and 
frequency. Fig, 10b typifies the case in which 
the shock 1$ near the trailing edge. Minimus 
values of Vf occur for this condition and over 
a small range of angle-of-attack there is little 
effect upon damping or flutter frequency. In 
fig. loc, the shock Is still further aft and, 
for the larger values of n, has moved off of the 
tralUng-edge. In this case, the effect of 
Increasing a;ig1e-of-attack reverses, with 
Increasing damping resulting for a > 0.5 deg. 
Also, the locus for a • 0.75 deg, indicates the 
cause of the multiple values of flutter 
boundaries shown in fig. 7. None of these 
mechanisms appear to address the experimental 
condition described by Yates et al33 and 
Houwlnk et al34 in which the flutter boundary 
was multiple valued for Mach numbers less than 
that at the bottom of the transonic dip. This 
observation led to the Investigation described 
In the following section of the effect of 
aeroelastic twist due to the static pitching 
moment, 


Effect of Stat ic pitching Moment - In (3), 
CntQ acFs’TrTTirel old wh1 ch “'is "adjusted to 
maintain the airfoil at a desired steady angle 
of attack. An alternative viewpoint Is adopted 


1n this section by Introducing the wind-off 
angle-of-attack, oq, and rewriting the static 
pitching moment equation as 


k^(n-ag) - l/2pu2(2b)2c^(c«, M) (11) 


which may be reorganized as 


V^« «r2 («.«^)/2 c,„(a, M) (12) 


The effects of angle-of-attack upon the 
flutter boundaries of the NACA 64A010 and MBB 
A-3 airfoils for Example 1 are shown In figures 
8 and 9. Changes in a of 1,5 degrees can Induce 
a 60S drop In Vf for the NACA 64A010 between 
Mach 0.775 and 0.80 and a 60< drop In Vf for 
the MBB A-3 between Mach 0.77 and 0.79. The 
reduced frequencies at flutter for the two 
examples range from kf » 0.12 for the higher 
Vf values to kf 0.2 for the lowest, A key 
feature is that the boundaries af fig. 8 show a 


Equation (121 Is a nonlinear equation relating 
the equilibrium angle-of-attack, a, to the speed 
Index, V, for given values of oq and M. With 
reference to strip-theory type analyses of 
wings, «o nwy be regarded as a "wing root 
angle-of-attack" and a as the local section 
angle-of-attack. The static pitching moment 
coefficient, (:„,(«), is plotted versus a and M 
for the NACA 64A010 and MBB A-3 airfoils in 
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figs, 11 «nd 12. For ». given Hich numNr. the 
pitching Mownt curves display three 
characteristics as a increases: 1)A range of a 
in which varies relatively linearly with a 
which corresponds to the transition from 
shockless flow to flow with mild shocks (see 
figs. S and 6), 2)a range of n 1n which Cm 
rises steeply corresponding to strong shocks 
moving aft on the airfoil, and 3)a range in 
which Cm again varies linearly with a with a 
slope miich is independent of H corresponding to, 
supersonic flow over the upper surface of the 
airfoil. Note particularly that the slope of 
the C||| curve in the first region varies 
gradually for the Sylnmetrical NACA 64A010 with 
both H and a whereas the slope of the Cm curve 
for the supercritical HBB A-3 is almost 
independent of H and a in this region. 

Figures 8 and 9 give flutter boundaries as 
a function of n. In order to determine flutter 
boundaries as a function of ao, taking into 
account twisting due to the steady pitching 
moment, the data from figs. 6 and 9 were cross.- 
plotted versus « for fixed M and solutions of 
(12) superimposed. Figures 13 and 14 give such 
plots for the NACA 64A010 at M • .8 and the MBB 
A-3 at M ■ .775. Intersections of the solution 
of (12) for a given oo with the flutter bound- 
ary curve represent flutter points at the indi- 
cated value of a. Figure 13 is typical for the 
NACA 64A010 in that only one flutter point 
occurs for each value of ao« in contrast, 
fig. 14 shows that the MBB A-3 has three flutter 
points for a range of bq near 4 deg. This 
occurs due to the steeper slope of the flutter 
boundary curves of the MBB A-3 coupled with the 
pitching moment behavior shown In fig, 12. 


Figures 15 and 16 present the effect of 
static pitching moment upon the flutter boundary 
of the two airfoils for 2° < oo < 6°. 

Comparing figs. 8 and 15, the effect for the 
NACA 64A010 is to steepen the flutter boundaries 
for Mach numbers between M » 0.75 and 0.8. 
Comparing Figs. 9 and 16, the effect for the MBB 
A-3 is much more pronounced. Multiple flutter 
points for a given Mach number cause the flutter 
boundary to curl back as V is decreased for oq 
> 3 deg. The amount of the curl -back in flutter 
Tiach number is similar to that shown in ref. 

33. At ao » 4 deg. Mf decreases 0.025 as 
Vf decreases from 1.1 to 0.65. The flutter 
boundaries shown in figs. 15 and 16 for Vf < 

0.8 correspond to Mach number and 
angle-of-attack combinations in which the shock 
has moved aft on the airfoil, in these cases, 
unmodeled boundary layer effects are probably 
Important, Nevertheless, the similarity of the 
flutter boundary curl back seen in fig, 16 to 
those shown in refs. 33 and 34 indicates that 
static aeroelastic twisting can have a 
significant effect upon transonic flutter. 

The multiple valued flutter curves shown in 
fig. 16 have an interesting interpretation in 
terms of flutter testing. At M-0.775 and 
Qo ■ 3 deg., the bendiiig mode root locus as a 
function of speed index 1s shoMi in fig. 17. 

The Speed index, V, and actual angle-of-attack, 
0 , are noted along the locus and show that the 


airfoil is twisted nose down as V is increased. 
Shown near V • 0,6 and a » 1 deg. is a local 
minimum in damping due to the proximity to the 
ao ■ 3 deg flutter boundary shown in fig. 16, 

As V increases, the mode becomes damped again 
and finally flutters at Vf * 1.3 where a • 

-0.7 deg, A small increase in either M or ao 
would lead to flutter at Vf » 0.6. The flutter 
point at u •> 133 rad/sec corrasponds to a 
classical type of flutter In which significant 
frequency coalescence occurs. The incipient 
flutter condition at u >> 85 rad/sec has much 
less frequency coalescence and is similar to 
cases which have been termed "single degree of 
freedom" or "shock Induced" flutter, 3“ 
Inspection of the static pressure distributions 
corresponding to these two conditions in fig, 6 
Indicates that the lower frequency flutter 
occurs when a shock has developed on the upper 
surface. 


Concluding Remarks 

The transonic small perturbation equation 
has been coupled with the structural equations 
of motion of a pitching and plunging airfoil and 
timennarchlng transient flutter solutions have 
been obtained. Accurate frequency and damping 
estimates were obtained by means of a complex 
exponential least squares curve fit of the re- 
sponses. Ine accuracy of the ttmewnarching 
calculations was established by comparison of 
results from the linearized transonic equation 
and by complrison of the flutter boundary ob- 
tained with the nonlinear equation with pub- 
lished results. Seven candidate numerical 
integration algorithms for the structural equa- 
tions were evaluated. The preferred algorithm 
is a modified state transition matrix integrator 
which was more accurate and stable for ’arger 
time steps than the others. 

The flutter boundaries of symmetrical NACA 
64A010 and supercritical MBB A-3 airfoils were 
determined for an example demonstrating a pro- 
nounced transonic dip. Comparison of flutter 
boundaries calculated using several nonlinear 
transonic aerodynamic codes show good agreement 
in predicting the transonic dip. The response 
at flutter was surprisingly linear in amplitude 
for angle-of-attack oscillations of up to three 
degrees. The effect of angle-of-attack upon the 
flutter boundaries of the two airfoils was 
determined. Changes of 1.5 deg in 
angle-of-attack can cause a 50( decrease In 
flutter speed index for the NACA 64A010 and a 
60% decrease for the MBB A-3. The slope of the 
flutter boundary with Mach number is larger for 
the MBB A-3 airfoil and appears to correlate 
with steady shock strength and locations on the 
airfoil. When aeroelastic twisting due to the 
static pitching moment is included, the steeper 
flutter boundary of the MBB A-3 leads to the 
occurrence of multiple flutter points for a 
given Mach number and a situation in whic'n the 
flutter boundary Mach number decreases as speed 
index is decreased. This effect has been 
observed in flutter model tests. 
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Fig. 1 Comparison of analytic and finite- 
difference calculations of lift 
coefficient due to pitching oscillation 
about midchord. 



TIME, «C 


Fig. 2 Fit Of transient response of MBB A-3 
airfoil at M X 0.8, a ■ 0 deg., and 
Vf » 0.645 for Example 1. 



MBB -A3 


Fig. 3 Airfoil profiles. 
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Fig. 4 Comparison of linearized root loci for 
Example 1 at M " 0.8. 
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Fig. 5 NACA 64A010 Steady pressure distri- 
butions as a function of Mach 
number and ang1e-of-altfck, 





Fig, 8 Effect of angle-of-attack upon NACA 

64AQ10 flutter boundaries obtained from 
HYTRAN2. 


Fig. 6 MBB A-3 steady pressure distributions as 
a function of Hach miiiber and angle-of- 
attack. 
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Fig. 9 Effect of «pgle-of-attack upon MOB A-3 
flutter bounUarles obtained from 
HVTRAN2. 




(a) (b) (c) 


F1g. 10 NACA 64A010 flutter mode root locus as 
a function of M, V, and 
angle-of-attack. 



Fig. 11 Static pitching moment of NACA 64A010 
^ibout x/b « -2 as a function of Mach 
number and angle-of-attack. 


Fig. 12 Static pitching moment of MBB A-3 about 
x/b « -2 as a function of Mach number 
and angle-of-attack. 



Fig. 13 Graphical determination of NACA 64A010 
flutter conditions at M » 0,8 Including 
effect of static pitching moment. 
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F1g, 14 Graphical determination of MBB A-3 
flutter conditions at M « 0,775 
Including the effect of static pitching 
moment. 




Fig, 15 The effect of static pitching moment on 
NACA 64A010 flutter boundaries. 


Fig, 17 MBB A-3 flutter mode root locus vs. V 
for Og » 3 deg. and M » 0.775. 



